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. Abstract. We describe a newly developed code for the extrapolation of nonlinear force- 

free coronal magnetic fields in spherical coordinates. The program uses measured vector 
, magnetograms on the solar photosphere as input and solves the force-free equations in 

the solar corona. The method is based on an optimization principle and the heritage of 
the newly developed code is a corresponding method in Cartesian geometry. We test the 
newly developed code with the help of a semi-analytic solution and rate the quality of 
our reconstruction qualitatively by magnetic field line plots and quantitatively with a 
' number of comparison metrics. We find that we can reconstruct the original test field with 

high accuracy. The method is fast if the computation is limited to low co-latitudes (say 
30° < 6 < 150°), but becomes significantly slower if the polar regions are included. 



> 

<N 

CN ' 1. Introduction 

\Q , 

, The solar corona is dominated by the coronal magnetic field because the 

magnetic pressure is several orders of magnitude higher than the plasma 
pressure. Knowledge regarding the coronal magnetic field is therefore im- 
portant to understand the structure of the coronal plasma and to get in- 
5_i ■ sights regarding dynamical processes like flares and coronal mass ejections. 

The direct measurement of magnetic fields is very difficult and such mea- 
surements are only occasional available, see, e.g., Lin et al. (2004). Well 
established are measurements of the photospheric magnetic field with the 
help of line-of-sight magnetographs (e.g., SOHO/MDI, NSO/Kitt Peak) or 
5_i ' vector magnetographs (e.g., currently with the Solar Flare Telescope/NAOJ, 

Imaging vector magnetograph/MEES Observatory, and in near future with 
SOLIS/NSO, Hinode, SDO/HMI). These photospheric fields can be extrap- 
olated into the solar corona by making suitable model assumptions. The 
simplest approach is to assume current-free potential fields, which can be 
computed from the photospheric line-of-sight magnetic field alone. Such 
source surface potential field models (Schatten et al., 1969; Neugebauer et 
al., 2002; Schrijver and Derosa, 2003) give already some insights about the 
coronal magnetic field structure, e.g., regarding the location of coronal holes 
and active regions. They cannot, however, be used to estimate the magnetic 
energy and helicity of the corona, in particular not the free energy available 
for eruptive phenomena. For reliable estimates of these quantities currents 
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have to be included. To lowest order the effects of plasma pressure and 
gravity can be neglected and we can assume that the currents are parallel 
to the magnetic field, the force-free assumption. A popular simplification of 
force-free fields are linear force-free fields (Chiu and Hilton, 1977; Seehafer, 
1978) where the electric current flow is parallel to the magnetic field with a 
global constant of proportionality a. This approach is in particular popular 
because linear force-free models require only line-of-sight magnetograms as 
input and contain only a single free parameter a, which can be specified 
by comparing magnetic field line plots with coronal images (Carcedo et al., 
2003; Marsch et al., 2004). In general a changes in space and taking this 
into account requires the use of nonlinear force-free models. A comparison of 
observationally inferred 3D magnetic field structures in a newly developed 
active regions (Solanki et al., 2003) with different extrapolated field models 
by Wiegelmann et al. (2005) revealed that linear force-free fields are better 
than potential fields, but nonlinear force-free models are even more accurate. 
The computation of nonlinear force-free fields is more challenging, how- 
ever for several reasons. Mathematically problems regarding the existence 
and uniqueness for various boundary value problems dealing with nonlinear 
force-free fields are still open, see Amari, Boulmezaoud, and Aly (2006) for 
details. Another issue is the numerical analysis for a given boundary value 
problem. An additional complication is to derived the required boundary 
data from photospheric vector magnetic field measurements. High noise in 
the transversal components of the measured field vector, ambiguities regard- 
ing the field direction and nonmagnetic forces in the photosphere complicate 
the task to derive suitable boundary conditions from measured data. 

Different approaches have been proposed for the nonlinear force-free ex- 
trapolation of vector magnetograms: 

— Upward integration method, (e.g., Wu et al., 1990, Cuperman, De- 
moulin, and Semel, 1991, Demoulin, Cuperman, and Semel, 1992, Amari 
et al., 1997). 

— Grad-Rubin-like method, (e.g., Grad and Rubin, 1958, Sakurai, 1981, 
Amari et al., 1997, Amari, Boulmezaoud, and Mikic, 1999, Wheatland, 
2004, Amari, Boulmezaoud, and Aly, 2006). 

— Different MHD relaxation methods, (e.g., Mikic and McClymont, 1994, 
Roumeliotis, 1996, Valori et al., 2005). 

— Green's function like method, (e.g., Yan and Sakurai, 2000, Yan and Li, 
2006). 

— Optimization method, (e.g., Wheatland et al., 2000, Wiegelmann, 2004). 

For a more complete review on existing methods for computing nonlinear 
force- free coronal magnetic fields see the review papers by Amari et al. (1997) 
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and Schrijver et al. (2006). The Grad- Rubin method as described in Amari 
et al. (1997) and Amari, Boulmezaoud, and Mikic (1999) has been applied to 
investigate particular active regions in Bleybel et al. (2002) and a comparison 
of the extrapolated field with 2D projections of plasma structures as seen in 
Ha, EUV and X-ray has been done in Regnier, Amari, and Kersale (2002) 
and Regnier and Amari (2004). The optimization code in the implementation 
of Wiegelmann (2004) has been applied to an active region in Wiegelmann et 
al. (2005a) and compared with Ha images. Wiegelmann and Inhester (2006) 
investigated the possibility to use magnetic field extrapolations to improve 
the stereoscopic 3D reconstruction from coronal images observed from two 
viewpoints. 

Recently Schrijver et al. (2006) compared the performance of six different 
Cartesian nonlinear force-free extrapolation codes in a blind algorithm test. 
All algorithms yield nonlinear force-free fields that agree well with the ref- 
erence field in the deep interior of the volume, where the field and electrical 
currents are strongest. The optimization approach successfully reproduced 
also the weak field regions and compute the magnetic energy content cor- 
rectly with an accuracy of 2%. In a coordinated study Amari, Boulmezaoud, 
and Aly (2006) obtained an accuracy of somewhat better than 2%. 

The good performance of the optimization method encourages us to de- 
velop a spherical version of the optimization code. The required full disk 
vectormagnetograms will become available soon (e.g., from SOLIS). The 
heritage of the newly developed code is a Cartesian force-free optimization 
method as implemented by Wiegelmann (2004). We outline the paper as 
follows. In Section 2 we describe our newly developed algorithm. Section 3 
contains a semi-analytic test case and the setup of computations to check 
the accuracy and performance of our code. We introduce figures of merit to 
rate the quality of our reconstruction in Section 4 and present the results 
of our test runs in Section 5. Finally we draw conclusions in Section 6 and 
give an outlook for future work. 



We solve Equations (1) and (2) with the help of an optimization principle 
as proposed by Wheatland et al. (2000) and generalized by Wiegelmann 
(2004). Until now the method has been implemented in Cartesian geometry. 



2. Method 



Force-free magnetic fields have to obey the equations 



(V x B) x B = 0, 
V B = 0. 



(1) 
(2) 
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Here we define a functional in spherical geometry: 

L= f \b~ 2 |(V x B) x B| 2 + [V • B| 2 j r 2 s'mO dr d9 dcj), (3) 
Jv 1 J 

It is obvious that the force-free Equations (1-2) are fulfilled when L equals 
zero. We normalize the magnetic field with the average radial magnetic field 
on the photosphere and the length scale with a solar radius. 

The functional (3) can be numerically minimized with the help of the 
iteration equations: 

where fi is a positive constant and 

f = v x (n a x b) - n a x (V x B) 

+V(f2 6 • B) - ft 6 (V ■B) + (ft 2 a + n 2 b )B (5) 

with 

ft a = B- 2 [(V x B) x B] (6) 

fl b = B- 2 [(V B) B] (7) 

The theoretical deviation of the iterative Equation (4) as outlined by Wheat- 
land et al. (2000) does not depend on the use of a specific coordinate system. 
Previous numerical implementations of this method have been done to our 
knowledge only in Cartesian geometry, however. Here we describe a newly 
developed implementation in spherical geometry. 

2.1. Implementation 

We use a spherical grid r, 6, 4> with n r , n$, grid points in radial direction, 
latitude 1 and longitude, respectively. Here we intend to compute the whole 
sphere r = lR s ... 2.57 R s , 6 = 0°... 180°, = 0°... 360°, but in principle 
one could limit the method also to parts of a sphere. To avoid the mathemat- 
ical singularities at the poles, we do not use grid points exactly at the south 
and north pole, but half a grid point apart 9 m i n = ^ and 6> max = 180° — 
The method works as follows: 

1. We compute the initial source surface potential field in the computa- 
tional domain from B r on the photosphere at r = lRg . 



1 9 corresponds to the co-latitude, with 9 — 0° and 9 — 180° at the south and north 
poles, respectively. 
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2. We replace Bg and at the bottom photospheric boundary at r = IRs 
by the measured vector magnetogram. 2 The outer radial boundary is 
unchanged from the initial potential field model. For the purpose of code 
testing we try additional other boundary conditions, see Section 3.2. 

3. We iterate for a force-free magnetic field in the computational box by 
minimizing the functional L of Equation (3) by applying Equation (4). 

4. The continuous form of Equation (4) ensures a monotonically decreasing 
functional L. For finite time steps, this is also ensured if the iteration 
time step dt is sufficiently small. If Lit + dt) > L(t) this step is rejected 
and we repeat this step with dt reduced by a factor of two. 

5. After each successful iteration step we increase dt by a factor of 1.01 to 
ensure a time step as large as possible within the stability criteria. This 
ensures an iteration time step close to its optimum. 

6. The iteration stops if dt becomes too small. As stopping criteria we use 
dt < 10~ 9 . 



3. Test Case 



2 It is as well possible to replace only parts of the photosphere (a limited region in 9 and 
4> direction) and to restrict the nonlinear force-free computation onto this region. This is 
in particular necessary if the observed photospheric vector magnetogram is only available 
for parts of the photosphere. In such cases the global magnetic field is basically a potential 
field and only locally (say in active regions) a nonlinear force-free field. For such limited 
regions in 9 and (j> we encounter the same problem as in Cartesian codes, that the lateral 
boundary conditions are unknown. One possibility is to describe the lateral boundaries 
with the help of a global potential field. The assumption of a potential field outside the 
computational domain restricts currents to the active region, but non current carrying 
field lines can leave the computational box. At the interface between the potential and 
non-potential field, a boundary layer as described in Wiegelmann (2004) for the Cartesian 
implementation of our code can be used. Full spherical force-free fields certainly do not 
have lateral boundaries and if one is interested in the details of interaction of two far apart 
active regions, it might be better to compute first a global low resolution force-free field 
(e.g., in future with SOLIS data) and then compute the field in the active regions with 
higher resolution (e.g., with Hinode data). Such an approach would allow also current 
carrying field lines to connect the two active regions, which might be important for the 
initiation of CMEs. 
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Figure 1, The synoptic map shows B r on the photosphere. We project the map onto the 
solar surface in Figure 2 with the disk center at Carrington Longitude 180°. 



3.1. Semi- Analytic Reference Field 

We test our newly developed code with the help of a known nonlinear force- 
free field model developed by Low and Lou (1990). The authors solved 
the Grad-Shafranov equation for axisymmetric force-free fields in spherical 
coordinates r, 9, (p. The magnetic field can be written in the form 



B 



1 



1 dA dA 



rs'mO \r 89 r dr 



e e + Qe, 



(8) 



where A is the flux function, and Q represents the ^-component of B, 
depending only on A. The flux function A satisfies the Grad-Shafranov 
equation 



d 2 A 1_ 

dr 2 



d 2 A 



dQ 



<V +Q dA 







where = cos 9. Low and Lou (1990) derive solutions for 



dQ 
dA 



a = constant 



by looking for separable solutions of the form 



(9) 



(10) 



11 



The solutions are axisymmetric in spherical coordinates with a point source 
at the origin. They have become a kind of standard test for nonlinear force- 
free extrapolation codes (Amari, Boulmezaoud, and Mikic, 1999; Wheatland 
et al., 2000; Wiegelmann and Neukirch, 2003; Yan and Li, 2006; Amari, 
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Figure 2. The figure shows the original reference field, a global potential field and the re- 
sults of a nonlinear force-free reconstruction with different boundary conditions (Case 1-4, 
see text). The color coding shows B r on the photosphere and the disk center corresponds 
to 180° longitude. 
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Figure 3. Evolution of L (as defined in 3) during the optimization process. The solid 
line corresponds to Case 1, the dotted line to Case 2, the dashed line to Case 3 and the 
dash-dotted line to Case 4. The horizontal dash-double-dotted line marks the discretization 
error of the original semi analytic solution. The top panel corresponds to the low resolution 
(20 x 40 x 80) and the bottom panel to the high resolution (40 x 80 x 160) computation. 
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Boulmezaoud, and Aly, 2006; Inhester and Wiegelmann, 2006; Schrijver et 
al., 2006) in Cartesian geometry because the symmetry in the solution is no 
longer obvious after a translation which places the point source outside the 
computational domain and a rotation of the symmetry axis with respect to 
the Cartesian coordinate axis. 

Here we use the Low and Lou solution in spherical coordinates. The 
original equilibrium is invariant in cj>, but we can produce a 3-D looking 
configuration by placing the origin of the solution with 1/4 solar radius 
offset to the sun center. The corresponding configuration is not symmetric 
in <j) anymore with respect to the solar surface as seen in the synoptic map 
in Figure 1 which shows B r on the photosphere. Let us remark that we use 
the solution for the purpose of testing our code only and the equilibrium is 
not assumed to be a realistic model for the global coronal magnetic field. We 
do the test runs on spherical grids (r, 9, <f>) of 20 x 40 x 80 and 40 x 80 x 160 
grid points. 

3.2. Boundary Conditions 

— Case 1: Boundary specified on the photosphere, source surface and in co- 
latitude at 9 = 30° and 9 = 150°. Optimization restricted in co-latitude 
to 30° < 9 < 150°. 

— Case 2: Boundary specified on photosphere, the source surface and in 
co-latitude at = f and 9 = 180 - f . 

— Case 3: Boundary specified on photosphere and on the source surface. 

— Case 4: Boundary specified only on the photosphere. This is the real- 
istic Case for real data, where measurements are only available on the 
photosphere. On the source surface the magnetic field is chosen from 
the initial potential field. 

For the computations done here we use a grid resolution of dO = d(f> = 4.5° for 
the low and d9 = d<f> = 2.25° for the high resolution test Case 3 . This means 
that the full spherical extrapolation cases (Case 2- Case 4) are limited in co- 
latitude by 2.25° < 9 < 177.75° and 1.125° < 9 < 178.875° for the low and 
high resolution computations, respectively. As the initial state we compute 
a source surface potential field in our computational domain. The source 
surface is a spherical shell where all field lines are assumed to become radial 
(Schatten et al., 1969). We locate the source surface at 1 + f ~ 2.57R S which 
is the outer radial boundary of our physical domain. The finite differences 

3 If the angle is given in radian we have even dr = d6 = d<f> for the test cases. The code 
allows, however, also the use of dr ^ dd ^ dcj>. 



sphere_ff_rev2.tex; 5/02/2008; 14:48; p. 9 



10 



Wiegelmann 



in <fi are cyclic. For Case 1 and 2 boundary values in 6 are specified and for 
Case 3 and 4 we interpolate the values at the poles. 

4. Figures of Merit 

In Table I we provide some quantitative measures to rate the quality of 
our reconstruction. Column 1 names the corresponding test case. Column 
2-4 show, how well the force and solenoidal condition are fulfilled, where 
column 1 contains the value of the functional L as defined in Equation (3) 
and L\ and L2 in column 3 and 4 correspond to the first (force-free) and 
second (solenoidal free) part of L. The evolution of the functional L during 
the optimization process is shown in Figure 3. Column 5 contains the Lqo 
norm of the divergence of the magnetic field 

|| V ■ B ||oo= sup |V ■ B| (12) 

and column 6 the norm of the Lorentz force of the magnetic field 

|| j x B ||oo= sup |j x B|. (13) 

The next five columns of Table I contain different measures which com- 
pare our reconstructed field with the semi-analytic reference field. These 
measures have been introduced by Schrijver et al. (2006) to compare a vector 
field b with a reference field B. 

— Column 7: Vector correlation: 

C vec = ^B i .b i /fc|B i | 2 ^|b i | 2 ] ' , (14) 

i \ i i / 

— Column 8: Cauchy-Schwarz inequality 

CCS = ^|£|' (15) 
where N is the number of vectors in the field. 

— Column 9: Normalized vector error 

S N = ^|b i -B i |/^|B i |, (16) 

i i 

— Column 10: Mean relative vector error 

^-^t (17) 
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Table I. 

The table rates the quality of our reconstructions with several figures of merit as explained 
in Section 4. We compute the figures for the whole sphere. For Case 1, where the 
computations have been limited in latitude, the reference field was specified in the cones 
6 < 30° and 6 > 150°. 



Model 


L 




L 2 


V ■ B Hex, 


II j XB ||oo 


Cvcc 


Ccs 


En 


Em 




Steps 


Time 




Spherical grid 20 x 40 x 80 




















Original 


0.04 


0.03 


0.01 


0.48 


1.37 


1 


1 








1 






Potential 


0.19 


0.006 


0.13 


4.61 


0.90 


0.66 


0.77 


0.71 


0.68 


0.75 






Case 1 


0.014 


0.010 


0.004 


0.42 


0.64 


0.9998 


0.9995 


0.012 


0.016 


1.008 


1889 


2 min 


Case 2 


0.013 


0.010 


0.003 


0.42 


0.64 


0.9998 


0.9992 


0.016 


0.023 


1.007 


31129 


23 min 


Case 3 


0.23 


0.20 


0.03 


3.00 


11.30 


0.998 


0.998 


0.04 


0.05 


0.99 


23824 


18 min 


Case 4 


0.48 


0.34 


0.14 


3.11 


11.30 


0.993 


0.95 


0.14 


0.29 


0.95 


10167 


7 min 




Spherical grid 40 x 80 x 160 




















Original 


0.005 


0.003 


0.002 


0.38 


0.71 


1 


1 








1 






Potential 


0.30 


0.0003 


0.30 


0.44 


0.23 


0.67 


0.77 


0.70 


0.67 


0.75 






Case 1 


0.0016 


0.0010 


0.0006 


0.38 


0.32 


0.99998 


0.99989 


0.004 


0.007 


1.002 


14522 


lh 26min 


Case 2 


0.0027 


0.0020 


0.0007 


0.40 


0.53 


0.99992 


0.9995 


0.011 


0.019 


0.9989 


672281 


65h 19min 


Case 3 


0.24 


0.20 


0.04 


6.63 


22.47 


0.996 


0.99 


0.086 


0.12 


0.96 


171143 


16h 57min 


Case 4 


0.39 


0.27 


0.12 


6.60 


22.47 


0.99 


0.93 


0.17 


0.32 


0.92 


120742 


12h OOmin 



— Column 11: Total magnetic energy of the reconstructed field normalized 
with the energy of the input field 

V- lb-1 2 

e - §HF- (18) 

The two vector fields agree perfectly if C vec , Ccs and e are unity and if £n 
and Em are zero. Column 12 contains the number of iteration steps until 
convergence and column 13 shows the computing time on 4 processors. 



5. Results 

5.1. Qualitative Comparison 

In Figure 2 we compare magnetic field line plots of the original model field 
(Original) with a corresponding potential field (Potential) and nonlinear 
force-free reconstructions with different boundary conditions (Case 1 - Case 
4). The colour coding shows the radial magnetic field on the photosphere, 
as also shown in the synoptic map in Figure 1. Carrington longitude 180° 
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corresponds to the disk center in Figure 2. The images show the results of 
the computation on the 20 x 40 x 80 grid. 

A comparison of the original reference field (Original in Figure 1) with our 
nonlinear force- free reconstructions (Cases 1-4) shows that the magnetic field 
line plots agree with the original for Case 1-3 within the plotting precision. 
Case 4 shows some deviations from the original, but the reconstructed field 
lines are much closer to the reference field than to the initial potential field 
(Potential). We cannot expect a perfect agreement with the reference field for 
Case 4, because here the outer radial boundary conditions are taken from the 
initial potential field model, which are different from the outer boundary of 
the reference field. It is well known from computations in Cartesian geometry 
that the lateral and top boundaries influence the solution in the box 4 . 

5.2. Quantitative Comparison 

The visual inspection of Figure 2 is supported by the quantitative criteria 
shown in Table I. For Cases 1 and 2 (where the boundary conditions have 
been specified on the photosphere, the source surface and in latitude) the 
formal force-free criteria (L, L\,L-i) are smaller than the discretization error 
of the analytic solution and the comparison metrics show an almost perfect 
agreement with the reference field. For Case 3 we still find a very good 
agreement between the reference field and our reconstruction and for Case 
4 the agreement is still reasonable. We are in particular able to compute 
the magnetic energy content of the coronal magnetic field approximately 
correct. A potential field reconstruction does obviously not agree with the 
reference field. The figures of merit show that the potential field is far away 
from the true solutions and contains only 75% of the magnetic energy. 

For a practical use of any numerical scheme the computing time certainly 
matters. As seen in the last column of Table I the computing time is quite fast 
if we restrict the computational domain to low co-latitudes (30° < < 150°, 
Case 1) but increases significantly for full sphere computations. The time 
needed for one iteration time step is approximately constant (0.04s on the 
low and 0.36s on the high resolution grid.), but the number of iteration steps 
needed until convergence increases by more than one order of magnitude if 
high co-latitude regions are included in the optimization. The reason is that 
the iteration time step dt, which adjusts automatically in our code, becomes 
much smaller and consequently the convergence speed becomes slow. This 
is well visible in Figure 3 which shows the evolution of the functional L 
with the number of iteration steps. For the low latitude computations (Case 
1, solid line in Figure 3) the functional L is much steeper than for the full 

4 In the Cartesian case this effect can be reduced by choosing a well isolated active 
region surrounded by a sufficiently large area with low magnetic field, see e.g., test Case 
II in Schrijver et al. (2006). 
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sphere computations (Cases 2-4, dotted,dashed, dash-dotted lines). The slow 
convergence of the full sun computations is directly related to the smaller 
time step and the time step itself is restricted by the physical grid resolution 
5 , which becomes very small close to the poles. We will address this point in 
the discussion (Section 6) and outline possible solutions to overcome these 
difficulties. 

Some insights regarding the performance of our newly developed code 
might be given by a comparison of the figures in Table I with the results of 
computations in Cartesian coordinates, as shown in Schrijver et al. (2006) 
Table I, where row (b) contains the results of our Cartesian optimization 
code. The upper part of Schrijver's Table I corresponds to a test Case 
where all six boundaries of a Cartesian box have been specified (Schrijver's 
Case I, which is somewhat similar to our Cases 1-3) and the lower part of 
Schrijver et al. (2006) Table I corresponds to Si CclSG where only the bottom 
boundary of a Cartesian box has been specified (Schrijver's Case II, which 
is equivalent to our Case 4 in spherical geometry). Both for Cartesian and 
spherical computations the correspondence with the original field is reduced 
if the boundary conditions are only specified on the photosphere. For real 
observed vector magnetograms we certainly have only photospheric data and 
it is therefore important to get a reasonable nonlinear force-free magnetic 
field reconstruction for this case. The errors of the different comparison 
metrics are still small and in a comparable range for Cartesian and spherical 
computations. The vector correlation is better than 99% and we got the 
magnetic energy correct within a few percent for all examples investigated 
in this paper. 



6. Discussion and Outlook 

Within this work we developed a code for the nonlinear force-free com- 
putation of coronal magnetic fields in spherical geometry. The method is 
based on an optimization principle. We tested the performance of the newly 
developed code with the help of a semi-analytic reference field. We find 
that the spherical optimization method works well and the accuracy of the 
reconstructed magnetic field configuration is comparable with the perfor- 
mance of a corresponding code in Cartesian geometry. The computation 
is reasonably fast if we limit the computation to low latitude regions, but 
becomes significantly slower if polar regions are included. Such a behaviour 
is well known for computations in spherical coordinates, see e.g. Kageyama 
and Sato (2004). The reason is that the physical grid converges to the poles. 
A fair approximation for the iteration time step is dt oc A 2 , where A is 

5 This condition is similar to the CFL-condition for time dependent problems. 
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the physical grid resolution. As the finest grid resolution restricts the time 
step, it becomes much smaller if high altitude regions (where A becomes 
very small) are included in the optimization. Our code has an automatic 
step size control and we find that dt automatically adjusts to much smaller 
values for full sphere computations. 

A possible solution is the use of the so called Yin- Yang grid, as developed 
by Kageyama and Sato (2004) and as used, e.g., for Earth mantle convection 
by Yoshida and Kageyama (2004). The Yin- Yang grid is composed of two 
identical complimentary grids, which partly overlap and together cover the 
solar surface with a quasi-uniform grid spacing (see Kageyama and Sato 
(2004) for details). Naturally the Yin- Yang grid has no poles at all and a 
considerably large iteration step can be used. There is certainly a numeri- 
cal overhead for transformations between the complimentary grids, but the 
limitations regarding the time step in polar regions for traditional spherical 
grids are more time demanding. For solar applications one has to consider, 
however, that photospheric magnetic field measurements are currently less 
accurate close to the poles. It might therefore be acceptable to restrict 
a nonlinear force-free computations onto equatorial regions, say between 
about 30° and 150° latitude 6 . For these low latitude regions the spherical 
optimization code as described here is reasonably fast. For application to 
observed vector magnetograms one has to consider, that the measured mag- 
netic field is not necessary force-free in the photosphere and in particular 
the transverse components of the magnetic field vector contain much noise. 
These problems are, however, also present for force-free extrapolations in 
Cartesian geometry and a preprocessing of the photospheric vector magne- 
tograms as described in Wiegelmann, Inhester and Sakurai (2006) helps to 
overcome these difficulties. 

The code solves the nonlinear force-free equations in the bounded domain 
between IRs and the source surface at 2.57Rs. The outer boundary is kept 
fixed from the initial potential field. All current carrying field lines have 
to close inside the volume. The domain outside 2.57Rs is not included in 
the model, because the force-free approach is not justified anymore here. 
A further step towards a consistent modelling of the solar corona would 
be the inclusion of non-magnetic forces, like plasma pressure, gravity and 
the dynamic pressure of the solar wind. This is in particular useful for 
several solar radii long structures like helmet streamer, where the plasma (3 
becomes finite (see e.g., Guo and Wu (1998) and Wiegelmann et al. (1998)). 
Neukirch (1995) found a special class of magnetohydrostatic solutions which 
are separable in spherical coordinates. Petrie and Neukirch (2000) extended 
the linear force-free Green's function methods towards the inclusion of non- 



6 Alternatively one could use a larger grid spacing for the computation of a force-free 
field in polar regions 
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magnetic forces in Cartesian geometry. Both approaches assume a kind of 
global linear force-free parameter a for the parallel part of the electric cur- 
rents. This is certainly a too restrictive condition to include the detailed 
information provided by measured vector magnetograms. The optimization 
method has been generalized for this aim in Wiegelmann and Inhester (2003) 
and implemented in Cartesian geometry. A corresponding implementation 
into spherical geometry is straightforward, but does certainly require addi- 
tional observational constraints, e.g., the tomographically reconstructed 3D 
coronal density distribution. 
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